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19.  Abstract 

beams  defined  from  vicinity  rays  will  exhibit  a  much  slower  breakdown  in  accuracy 
as  the  scale  length  of  the  medium  given  by  v  approaches  the  beamwidth. 

Vv 

Since  second  spatial  derivatives  of  velocity  arc  not  required  by  the  new  technique, 
parameterization  of  the  medium  is  simplified,  and  reflection  and  transmission  of 
beams  can  be  calculated  by  applying  Snell's  law  to  both  vicinity  rays  and  central 
rays . 
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TECHNICAL  SUMMARY 


The  objective  of  this  project  is  to  determine  the  yield  bias  of  underground  nu¬ 
clear  tests  induced  by  the  presence  of  a  high  velocity  descending  slab  beneath 
the  test  site.  Specifically,  the  effect  of  the  Aleutian  slab  is  being  investigated 
on  the  US  underground  tests  Longshot,  Milrow,  and  Cannikan.  P  wave  seis¬ 
mograms  will  be  synthesized  using  dynamic  ray  tracing  and  superposition 
of  Gaussian  beams  in  three-dimensional  models  of  the  Aleutian  slab  deter¬ 
mined  from  P  travel  time  delays.  Focusing  and  defocusing  and  multipathing 
at  teleseismic  distances  will  be  evaluated  by  comparison  of  observed  with 
synthetic  seismograms  of  the  Aleutian  tests. 

The  slab  problem  requires  that  accurate  body  waves  be  synthesized  for 
three  dimensional  source- receiver  geometries  in  a  three-dimensionally  vary¬ 
ing  Earth  model.  The  rapid  variation  of  Earth  structure  in  the  vicinity  of 
the  slab  requires  care  in  assessing  the  validity  and  appropriateness  of  apply¬ 
ing  asymptotically  approximate  techniques  of  synthesizing  body  waveforms. 
During  the  first  six  months  of  the  project,  we  have  developed  an  improved 
method  of  dynamic  ray  tracing  and  Gaussian  beam  synthesis  for  application 
to  this  problem.  The  improved  method  does  not  make  any  paraxial  approxi¬ 
mations  and  defines  bearnwidth  in  terms  of  a  Fresnel  volume.  By  avoiding  the 
paraxial  approximation,  we  seek  to  develop  a  beam  method  that  will  remain 
valid  in  rapidly  varying  Earth  models  up  to  the  point  at  which  asymptotic 
ray  theory  fails.  The  conventional  dynamic  ray  tracing  and  beam  superposi¬ 
tion  method  breaks  down  as  the  beam  width  approaches  the  scale  length  of 
the  medium,  which  may  be  long  before  ray  theory  itself  breaks  down. 

The  Fresnel  bearnwidth  criterion  has  been  tested  for  diffraction  from  caus¬ 
tics  by  comparison  with  the  predictions  from  WKBJ  seismograms.  Although 
the  diffraction  predicted  by  WKBJ  seismograms  is  also  an  approximation,  it 
is  appropriate  to  compare  the  two  techniques  since  both  techniques  are  devel¬ 
oped  using  an  equivalent  level  of  asymptotic  approximation.  Good  agreement 
is  obtained  between  the  modified  beam  method  and  WKBJ  seismograms  for 
a  simple  model  consisting  of  a  gradient  discontinuity. 
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1  Introduction 


Many  high  frequency  asymptotic  solutions  of  the  wave  equation  have  been 
developed  as  effective  tools  for  computing  wave  fields  in  inhomogeneous  three 
dimensional  media.  Two  of  the  most  widely  applied  are  the  WKBJ/Maslov 
method  (Chapman,  1978;  Chapman  and  Drummond,  1982)  and  the  Gaus¬ 
sian  beam  method  (Babich  1980;  Popov,  1982;  Cerveny  et  al.,  1982;  and 
Cerveny  and  Psencik,  1983).  Both  of  these  techniques  estimate  the  kine¬ 
matic  and  dynamic  properties  of  a  wavefront  from  approximate  solutions 
to  the  elastodynamic  wave  equation  based  on  ray  theory.  The  superposi¬ 
tion  techniques  of  Gaussian  beams  and  YVKBJ  plane  waves  as  well  as  their 
stationary  phase  approximation  in  geometric  ray  theory  al!  require  similar 
amplitude  and  weighting  functions.  These  amplitude  functions  can  be  found 
by  integrating  a  system  of  equations  known  as  the  dynamic  ray  tracing  equa¬ 
tions. 

The  dynamic  ray  tracing  equations  can  be  derived  from  either  the  eikonal 
equation  by  substitution  of  a  paraxial  approximation  (Cerveny  and  Hron, 
1980;  Cerveny,  1985),  or  from  the  parabolic  wave  equation  (Cerveny  et  ah, 
1982;  Popov,  1982;  Cerveny  and  Psencik,  1983).  The  dynamic  ray  tracing 
equations,  however,  have  both  limitations  and  complications.  The  limita¬ 
tions  are  associated  with  the  use  of  the  paraxial  approximation,  and  the 
complications  are  due  to  the  use  of  multiple  coordinate  systems. 

The  limitations  associated  with  the  paraxial  approximation  are  exhibited 
whenever  the  dynamic  ray  tracing  equations  are  used  to  estimate  the  travel 
time  and  amplitude  at  a  point  in  the  neighborhood  of  ray  from  a  second  order 
Taylor  expansion  of  the  wavefront  at  a  point  along  the  ray.  The  Taylor  ex¬ 
pansion  is  the  essential  step  in  the  definition  of  Gaussian  beams  and  paraxial 
rays.  The  region  in  which  the  error  of  this  Taylor  expansion  remains  below 
some  specified  threshold  is  generally  referred  to  as  the  paraxial  vicinity.  The 
fundamental  problem  with  the  paraxial  approximation  is  that  it  is  not  simple 
to  specify  the  spatial  bounds  of  the  paraxial  vicinity  in  a  three-dimensionally 
varying  model.  In  general,  one  must  not  attempt  to  evaluate  the  Taylor  ex¬ 
pansion  too  far  from  the  central  ray,  but  it  is  unknown  how  the  error  grows 
away  from  the  central  ray. 

The  complications  associated  with  the  use  of  two  coordinate  systems  are 
best  appreciated  by  considering  the  most  general  case  of  a  three-dimensionally 
varying  medium.  In  three-dimensionally  varying  media,  the  usual  approach 


is  to  specify  the  dynamic  ray  tracing  equations  using  two  coordinate  systems: 
ray  coordinates,  usually  consisting  of  the  take-off  angles  at  the  source,  and 
ray  centered  coordinates,  consisting  of  a  orthogonal  curvilinear  system  that 
moves  along  with  the  ray  (figure  1).  The  use  of  two  coordinate  systems, 
while  having  the  advantage  of  converting  a  non-linear  Ricatti  equation  into 
a  system  of  linear  equations,  increases  the  number  of  equations  needed  to 
describe  the  quantities  affecting  the  amplitude  of  the  wavefield.  In  either 
the  fixed  Cartesian  or  ray  centered  coordinates,  the  standard  dynamic  ray 
tracing  equations  require  the  specification  of  the  second  spatial  derivatives 
of  velocity  along  a  ray.  This  either  forces  the  model  to  be  parameterized 
with  continuous  first  derivatives  of  velocity  or  complicates  the  integration  by 
requiring  jump  conditions  on  the  dynamic  quantities.  These  jump  conditions 
are  obtained  by  matching  the  paraxially  estimated  phase  on  either  side  of 
the  discontinuity  in  gradient  (Cerveny  and  Ilron,  1980;  Cerveny,  1985). 

In  this  paper,  we  develop  a  new  system  of  dynamic  tracing  equations 
and  Gaussian  beams  that  eliminates  many  of  these  problems.  This  system 
is  derived  from  the  Hamiltonian  of  a  ray.  We  will  define  a  ’’vicinity  ray” 
to  be  a  ray  in  the  neighbourhood  of  the  central  ray.  Each  vicinity  ray  has 
slightly  different  initial  take-off  angles  with  respect  to  the  central  ray  (figure 
2).  The  locus  of  the  vicinity  ray  is  not  paraxially  estimated  from  the  standard 
dynamic  ray  tracing  equations,  but  rather  is  determined  much  more  precisely 
by  integrating  a  neiy  set  of  differential  equations,  which  we  refer  to  as  the 
’’vicinity  ray  tracing  equations.”  Gaussian  beams  are  defined  by  assigning  a 
Gaussian  distribution  of  amplitude  to  each  centra!  ray.  The  width  of  the 
Gaussian  is  taken  to  be  the  Fresnel  volume  surrounding  the  central  ray. 
Since  beamwidths  are  related  to  the  Fresnel  volume,  diffracted  wavefields 
can  be  accurately  estimated  by  a  superposition  of  Gaussian  beams  without 
the  ambiguity  associated  with  a  freely  varying  beamwidth  parameter. 

In  the  following  sections  we  first  review  the  derivation  of  the  standard 
dynamic  ray  tracing  system  and  the  limitations  of  the  paraxial  approxima¬ 
tion.  Next  the  vicinity  ray  tracing  system  is  derived  from  the  Hamiltonian 
of  a  ray,  in  which  no  paraxial  approximations  are  made.  Expressions  for  the 
travel  time  and  wavefront  curvature  in  the  neighborhood  of  a  central  ray  are 
derived  using  this  system.  Gaussian  beams  are  defined  using  vicinity  rays  to 
approximate  the  Fresnel  volume.  Finally,  seismograms  are  synthesized  and 
compared  in  a  simple  one-dimensional  model  using  the  WKBJ  method  and 
superposition  of  Gaussian  beams  defined  from  vicinity  rays. 
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2  Physical  and  mathematical  system 

Consider  an  arbitrary  ray  corresponding  to  a  P-wave  and  introduce  ray  cen¬ 
tered  coordinates  s,q\,qz  (figure  1).  The  orthogonal  ray  centered  coordinate 
system  along  the  central  ray  Q  and  its  computation  are  described  in  Popov 
and  Psencik  (1976),  Psencik  (1979),  and  Cerveny  and  Hron  (1980).  The 
ray  centered  coordinates  are  limited  to  a  vicinity  of  the  origin  ( qt  =  0)  in 
which  the  central  ray  field  is  regular.  In  figure  1,  the  coordinate  s  measures 
the  arclength  along  a  central  ray  from  an  arbitrary  reference  point.  ql  and 
q2  represent  length  coordinates  and  form  a  two-dimensional  Cartesian  coor¬ 
dinate  system  in  the  plane  normal  to  0  at  O,  with  origin  at  SI.  All  three 
components  (s,q^q2)  in  the  ray  centered  coordinate  system  depend  on  the 
azimuth  and  vertical  take-off  angle  (<t>,6).  The  basis  of  the  coordinate  system 
forms  a  right  handed  system  of  the  three  unit  vectors  t,  ey  and  e2  where  t  is 
the  unit  tangent  vector  to  the  central  ray  Q. 


2.1  Limitations  in  the  dynamic  ray  tracing  system 

2.1.1  The  paraxial  approximation 

The  standard  dynamic  ray  tracing  system  can  be  derived  from  either  the 
eikouai  equation  (Cerveny  and  llrou,  1980;  Madariaga,  1984;  Cerveny,  1985) 
or  from  the  parabolic  wave  equation  (Popov,  1982;  Cerveny  and  Psencik, 
1983).  In  either  derivation,  a  paraxial  approximation  is  assumed  at  some 
stage,  which  involves  a  Taylor  expansion  of  the  wavefield  about  the  central 
ray.  This  approximation  and  other  approximations  occuring  in  the  derivation 
are  not  always  applied  consistently  and  terms  arc  oimttc  d  without  specifying 
validity  conditions. 

To  illustrate  the  problems  with  the  dynamic  ray  tracing  system,  let  us 
review  the  derivation  of  the  two-dimensional  dynamic  ray  tracing  equations 
starting  from  the  eikonal  equation.  The  eikonal  equation  in  2-D  is 


h^d.*  wy  v2 


where  V  -  v(.v,q).  h  is  a  scale  factor  in  the  ray  centered  coordinates  and  will 
be  discussed  subsequently.  The  travel  time  of  a  vicinity  ray  r(.s,g)  can  be 
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approximated  at  q-0  (Cerveny  and  Hron,  1980;  Cerveny  and  Psencik,  1983; 
Cerveny,  1985)  by: 

r(.s,9)  %  t(s)  1  X~M(s)q 2  (2) 

where  —  0  and  M  - 

From  equation  (2),  it  follows  that 


dr{s,q)  _  dr  (a)  1  0M{s)  ,  1  1 

(Is  (9s  2  da  ^  v  2 


where  v  -■  v(s,0). 
order  terms  gives 


dr(9,q) 

Oq 


Mq 


(3) 


Substituting  (3)  into  equation  (1)  and  neglecting  higher 

,  V.,’ha/v  L  w 


Rearrangng  terms  in  (4)  gives 


M'W 


1 


V'2  v2/i2 


(5) 


By  expanding  the  velocity  V  to  second  order  terms  with  respect  to  q-0, 


V’  -  v(s,q)  ~  v  I  v,qq  I  -v^qq2 , 

the  right  side  of  equation  (5)  can  be  approximated  by 
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V2  v2h2 


(6) 


(7) 


(Cerveny  and  Hron,  1980;  Cerveny  and  Psencik,  1983;  Cerveny,  1985),  where 

d*  v 

v  —  —  r 
djf- 

The  standard  dynamic  ray  tracing  system  is  obtained  from  equation  (5)  by 
using  equation  (7)  and  assuming  h  i.  This  gives 

dM 


da 


4  vM2  i  ^ q 2 


0 


(8) 
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Since  the  derivation  of  the  dynamic  ray  tracing  system  includes  second 
order  terms,  any  omitted  terms  must  he  carefully  evaluated.  Consider  the 
scale  factor  h.  The  scale  factor  h  is  given  hy 

h  I  I  ---</  (9) 

V 

Because  it  is  assumed  h  1  in  equation  (5),  the  neglected  term  2'^  q  of  h2  in 
equation  (.r>)  must  be  vanishly  small,  i.e., 


The  condition  in  equation  (10)  describes  the  applicability  of  the  dynamic 
ray  tracing  system.  It  says  that  extrapolation  of  the  wavefield  away  from 
a  central  ray  using  the  paraxial  approximation  will  break  down  rapidly  as 
the  scale  length  of  the  medium  increases.  The  extrapolation  distance  must 
be  much  less  than  the  scale  length  of  the  medium.  For  Gaussian  beams,  it 
implies  that  the  beam  width  must  be  much  less  than  the  scale  length  of  the 
medium.  This  can  be  a  severe  restriction  in  rapidly  varying  models,  in  which 
the  criterion  for  validity  of  ray  theory  (  wavelength  <C  scale  length)  is  still 
well  satisfied. 

The  term  given  by  the  left  side  of  the  inequality  (10)  would  also  be  omit¬ 
ted  if  the  derivation  of  the  dynamic  ray  tracing  equation  had  instead  started 
from  the  parabolic  wave  equation.  In  this  case,  the  omission  of  this  term  can 
occur  through  the  lack  of  internal  consistency  in  deriving  this  system  from 
the  parabolic  wave  equation  or  from  the  cikonal  equation  (i.e.,  approxima¬ 
tions  of  h  that  are  inconsistent  with  estimates  of  asymptotic  order  given  by 
wavelength/scale  length). 


2.1.2  The  P  and  Q  matrices  and  coordinate  systems 


Equation  (8)  is  a  non-linear  ordinary  differential  equation  of  the  first  order 
Ricatti  type.  This  equation  can  be  solved  by  elementary  analytical  methods. 
Following  G'crveny  and  llron  (1980),  the  2  I)  system  given  by  (8)  can  be 
generalized  to  a  3  1)  system  of  linear  differential  equations  by  introducing  a 
2x2  matrix  M: 


A  / 


dQ 


't,q 


on 
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where  Q  is  a  2x2  matrix.  Define  a  2x2  matrix  1*  as: 


P  -  v 


i  dQ 
ds 


(12) 


By  substituting  equations  (11)  and  (12)  into  equation  (5),  the  dynamic  ray 
tracing  equations  in  3-1)  can  be  written  as 


dQ 

ds 


vP 


dP 

ds 


I 


SQ 


where  )  pt-  —  and  7  is  ray  parameter.  S  is  given  as: 


(13) 


The  dynamic  ray  tracing  system  has  8  equations  for  real  Q  and  P,  and 
16  for  complex  Q  and  P  in  3-D  and  is  specified  in  ray  centered  coordinates 
(si9i)92)  and  ray  coordinates  (71,72).  Cerveny  (1985)  has  shown  that  only  8 
equations  are  generally  needed  for  Gaussian  beams.  The  number  of  equations 
can  be  reduced  further  if  only  one  coordinate  system  could  be  used. 

The  standard  dynamic  ray  tracing  system  generally  will  have  oir  diagonal 
terms  in  the  matrices  Q  and  P.  The  existence  of  these  off  diagonal  terms  is 
due  to  the  use  of  tw-o  coordinate  systems  in  describing  the  equations. 


2.1.3  Beamwidths 

The  idea  of  beamwidth  is  somewhat  nebulous  in  standard  Gaussian  beam 
theory,  and  a  proper  mathematical  or  physical  meaning  of  complex  param¬ 
eters  in  the  Q  and  P  matrices  is  not  considered  in  routine  applications  of 
the  method.  Complex  Q  and  P  can  be  shown  to  be  a  consequence  of  an  ap¬ 
proximate  solution  for  complex  rays  emanating  from  a  source  having  a  small 
imaginary  part  to  its  location  in  space  (Kelscn,  198-1;  Wu,  1985).  In  practice, 
beamwidths  are  defined  somewhat  arbitrarily  and  are  adjusted  to  minimize 
errors  in  the  beam  superposition  (K limes,  1988;  Kim  and  Garmany,  1985)  or 
tuned  to  minimize  errors  associated  with  rapid  variations  in  velocity  (Weber, 
1988).  White  et  al.  (1987)  have  shown  that  optimum  beamwidths  strongly 


depend  on  the  specific  wave  propagation  problem  and  the  particular  type 
of  boundary  interactions  occuritig  in  the  problem.  One  of  reasons  why  the 
concept  of  optimum  beamwidths  does  not  work  well  is  that  the  total  energy 
of  each  beam  differs  for  different  initial  beamwidths.  This  is  true  for  all  of 
the  various  optimal  beamwidths  that  have  been  proposed.  If  energy  flux 
is  to  be  conserved  within  a  ray  tube,  then  a  normalization  condition  must 
be  applied  with  respect  to  the  different  initial  beamwidths.  The  following 
section  shows  how  these  many  of  these  problems  in  the  standard  Gaussian 
beam  method  can  be  remedied  by  using  a  dynamic  ray  tracing  system  de¬ 
rived  from  the  Hamiltonian  of  a  ray  and  applying  a  normalization  condition 
and  conservation  law  of  energy  flux  along  a  wavefront. 

2.2  Vicinity  ray  tracing  system 

2.2.1  Derivation 

Let  us  consider  the  high  frequency  asymptotic  solution  to  the  wave  equation 
in  an  inhomogeneous  medium.  In  order  to  obtain  the  desired  approximation, 
let  us  assume  that  the  displacement  u  is  expressed  in  the  following  form  in  a 
generalized  coordinate  system. 

w(qt ,a;)  A(q,)e (14) 

where  i  —  1,2, ...«  is  an  n-dimensional  configuration  space  whose  coordinates 
are  the  n  generalized  coordinates  </,.  A  and  r  are  an  amplitude  function  and 
a  phase  function.  Roth  A  and  r  are  functions  which  can  be  assumed  to  be 
slowly  varying  with  respect  to  the  wavelength  A.  The  Hamiltonian  of  a  ray 
is  applied  in  this  study  to  determine  the  functions  A  and  r. 

Since  Fermat’s  principal  of  least  time  can  be  expressed  by  the  operations 
of  variational  calculus,  a  Lagrangian  and  Hamiltonian  of  a  seismic  ray  can  be 
defined  similar  to  those  used  in  describing  the  mechanics  of  particles.  Kara 
and  Madariaga  (1988)  ,  for  example,  used  the  Hamiltonian  of  a  seismic  ray 
to  develop  a  perturbation  theory  to  compute  the  amplitude  and  travel  time 
of  a  vicinity  ray  with  respect  to  a  reference  ray. 

By  applying  Fermat’s  principle,  the  travel  time  r  from  source  (s0)  to 
receiver  (sr)  can  be  written  as  a  path  integral  over  the  Langrangian,  L,  by 

t  J'  L(q„q,)ds  (15) 
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where  <7«  =  The  Lagrangian  of  the  ray  in  generalized  coordinates  is  given 
by 

£(9i>92>9i>92>-’)  T  ‘(X^^9,2)‘/2  ('6) 

i  i 

where  A;  is  a  scale  factor  of  the  i-th  unit  vector  and  V  is  the  correspond¬ 
ing  velocity.  Ray  tracing  in  a  general  coordinate  system  (e.g.,  Cartesian, 
spherical,  and  cylindrical)  is  given  by  the  Euler  equation, 

d  dlj  dlj  . 

da  kidq,  A,<9<7,  ’ 

which  is  derived  from  the  variational  principal.  The  generalized  momentum 
Pi  associated  with  the  coordinate  qi  is  given  as  follows  (Goldstein,  1980, 
p.339). 

Pi  ^  (18) 

hldql 

The  terms  canonical  momentum  or  conjugate  momentum  are  also  used  for  p;. 
The  Hamiltonian  of  the  ray  in  the  ray  centered  coordinates  can  be  obtained 
from  the  Lagrangian.  The  Lagrangian  of  a  ray  in  ray  centered  coordinates  is 
obtained  from  equation  (16)  as 


T(?i,<72,<7i,92,s)  -  -jjl**2  t  q]  \  q\)XI 7 


where 


h  -  hx  =  1  +  Y.  f  9* 


and  where  qi  ■-  u,  -  and  ht  /13  1.  V  -  f(s, 91,92)  is  the 

velocity  of  a  vicinity  ray  and  v  -  t)(s,0,0)  is  the  velocity  of  a  central  ray. 
The  Lagrangian  in  equation  (19)  has  s  as  an  independent  variable.  The 
conjugate  momentum  p,  can  be  expressed,  by  using  equation  (18), 

dL  1  qt 

P'  ~  H  v  f  <i| 

Pl  ,  9L  (21) 

^2  '  Jh2  1  gf  \  q\ 
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Equation  (21)  can  be  solved  for  <ji  and  q2,  yielding 

_ hpi _ 

Vy/l  V2(pi2  <  P22) 

hp2 


0 1 


02 


V'v/l  -  v/2(p.2  t  P22)’ 

The  Hamiltonian,  H,  is  expressed  as  follows: 

//(<7i,<72,PhP2,s)  =  Pi <?i  1  P2<h  ‘  LiOxyOi'Oi'Oi'3)- 


(22) 


(23) 


By  substituting  equation  (22)  into  equation  (23),  the  Hamiltonian  of  a  ray 
in  the  ray  centered  coordinates  is  obtained, 


~  y\l  V!(pJ  i  pi)]'11- 


(24) 


The  vicinity  ray  tracing  system  in  the  ray  centered  coordinates  can  be 
described  in  terms  of  the  canonical  equations  from  the  Hamiltonian  defined 
in  (24).  The  canonical  transform  in  the  ray  centered  coordinates  is  obtained 
from  equation  (23). 

dqi  OH  h  V  pi 
da  0p\ 
dq2  <9// 
dp2 


A 

h  V  pi 


da 


(25) 


a  1 

S-|2 

II 

OH 

0q\ 

V 

^'9i  h  A 

V2 

h  V'*,  (Pi  f  P2) 
A 

dpi 

OH 

A  - 

^  ,4 

h  Vm  ( P ?  I  p\) 

da 

dqi 

-  V  A 

V2 

A 

(26) 


where  A  ~  \J\  -  V2(p\  +  p22).  Equations  (25)  and  (26)  are  comparable  to 
the  dynamic  ray  tracing  equations  (13),  but  no  paraxial  approximations  have 
been  made.  Cerveny  (1988)  has  also  briefly  described  the  derivation  of  equa¬ 
tions  (25)  and  (26),  showing  how  they  become  the  standard  dynamic  ray 
tracing  equations  if  paraxial  approximations  are  substituted  for  h  and  V  . 

To  simplify  equations  (25)  and  (26),  let  us  define  q  as  the  angle  difference 
between  the  tangential  vectors  of  a  central  ray  and  a  vicinity  ray  in  the  (f,ej) 
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plane  and  £  as  the  angle  difference  between  the  tangential  vectors  of  a  central 
ray  and  of  a  vicinity  ray  in  the  (t,e2)  plane  in  the  ray  centered  coordinates 
(figure  2).  qi  is  the  distance  from  the  central  ray  to  the  the  vicinity  ray  along 
the  ii.  From  figure  3  it  is  seen  that  the  curvature  (Ki)  of  the  wavefront  is  a 
function  of  tan  tj,  and  qi , 


tan  qt 


1 


(27) 


where  i  —  1,2  and  Rt  is  the  radius  of  curvature  of  the  wavefront.  Using  the 
definition  of  tj  from  figures  2  and  3,  equation  (25)  can  be  rewritten  as 


dqx 

da 

dq2 

da 


h  tan  tj 
h  tan  ( 


(28) 


Using  equations  (21)  and  (28),  pi  can  be  expressed  with  respect  to  tan  q  and 
tan  (. 

tan  rj 
Pi 


Pi 


VB 
tan  £ 

T7T 


(29) 


where  B  =  v^l "+  tan2q  +  tan2(, 

Differentiating  equation  (29  )  with  respect  to  s  yields 

dpi  aec2i]  2  d7/  \\,tan  tj  tan  ij  tan  £  svc2(d( 

~d7  =  VlP{[  '  tan  °ds  "W"  - VB* - ds 


sec2(,  ,  Y\,lan  £  tan  q  tan  (  scr2q  dij 

,lan’l)d. 


£  =v&K‘  vm  vip  i. 

The  quantity  A  can  be  rewritten  using  expressions  for  p,  given  in  (29): 


(30) 


A-\j\  v2(p]jrp\)  yi 


tan2rj  I  tan2£  1 

B2  “  Tl 


tj  and  (  are  located  within  t  cx  plane  and  within  t  e2  plane  respectively, 
and  both  planes  arc  orthogonal  to  each  other.  This  orthogonality  simplifies 
the  derivation  of  vicinity  ray  tracing  system,  jj  and  can  be  obtained  by 
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equating  equations  (26)  and  (30).  Finally,  the  vicinity  ray  tracing  system  in 
the  ray  centered  coordinates  can  he  given  as 

dqi  u  , 

.  -  ~  h  tan  j] 

as 


dq2 

ds 


h  tan  £ 


(31) 


-  =  cos2ij(--~B2tan  i)  I - — -  f  I)  tan  n  tan  C) 

ds  V  cos2tj 


d( 

ds 


V 

—  c.os2(,  (--*-  B2tan  (  |  --  )  C  tan  q  tan  £) 


D 

cos2( 


(32) 


where 


C  h...  ^ 


n  -  h  l-  B2 

u  —  n 

Note  that  C  and  D  and  the  equations  for  i)  and  £  depend  on  the  velocity  of 
the  medium  along  the  vicinity  ray,  V'(.«i,0,92)  or  F(.s,0,ft)  rather  than  on  the 
velocity  of  the  medium  along  the  central  ray,  ^($,0,0)  -  v(s).  For  a  velocity 
model  specified  in  Cartesian  coordinates,  the  velocity  V  and  its  derivatives 
V,,,  can  be  calculated  by  transforming  the  positions  of  the  vicinity  ray  in  ray 
centered  coordinates  (0,ft,0)  and  (0,0, ft)  to  Cartesian  coordinates. 

C  and  D  can  be  expanded  at  ft  —  0,  by  using  equations  (7),  (8)  and  (29). 


*»!<?! 


«7i 


•  i<n 


+v 


?1 


(tan2 t]  -f  tan2() 


D  %  q2  V«H±”™-Sl{tan2Tj  +  tan2 Q  (33) 

v  v 

There  is  no  advantage,  however,  in  making  such  an  expansion.  The  accuracy 
of  this  expansion  decreases  as  the  distances  q  increase,  and  it  requires  the 
calculation  of  the  second  spatial  derivative  of  velocity. 

Since  the  vicinity  ray  tracing  system  calculates  ft  and  ft  values  by  using 
V  and  V,qi  it  is  not  necessary  to  employ  the  method  of  matching  paraxial 
phase  (Cerveny  and  llron,  1980)  to  determine,  new  initial  conditions  on  ft 
and  ft  when  vicinity  rays  are  transmitted  through  or  reflected  by  discontinu¬ 
ities.  Since  second  spatial  derivatives  of  velocity  are  not  used  in  the  vicinity 
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ray  tracing  system,  no  jumps  <7,  or  1 7,  are  induced  by  velocity  gradient  dis¬ 
continuities.  At  first  order  discontinuities  in  velocity,  new  initial  conditions 
on  and  77,  are  computed  by  simply  applying  Snell’s  law  to  both  the  central 
ray  and  the  vicinity  ray. 

2.2.2  Initial  conditions 

The  initial  conditions  of  the  vicinity  ray  tracing  system  at  the  source  point 
depend  on  the  type  of  source.  For  a  point  source  the  initial  conditions  are 

Vi  I.  0  0 

Vi  I.-  0  v!  (34) 

and  for  a  line  source  they  are 

<k  l«  0  Vi 

Vi  l.-o"  Vi  (35) 

where  superscript  1  denotes  the  initial  value  of  the  parameter.  The  initial 
condition  on  qj  in  the  case  of  a  line  source  depends  on  the  intensity  or  shape 
of  the  source.  q(  will  be  the  half  length  of  the  line  source. 

When  the  wavefield  is  computed  by  a  superposition  of  Gaussian  beams, 
the  initial  value  chosen  for  will  depend  on  the  density  of  beams  in  the 
superposition.  Beamwidths  are  taken  to  be  roughly  equal  to  the  Fresnel  vol¬ 
ume  surrounding  the  central  ray.  The  Fresnel  volume  is  estimated  from  the 
frequency  and  the  separation  of  the  vicinity  rays  from  the  central  ray.  To 
achieve  an  accurate  estimate  of  the  Fresnel  volume,  the  spacing  of  vicinity 
rays  is  taken  to  overlap  the  spacing  of  central  rays.  In  the  synthetic  seismo¬ 
grams  shown  in  a  later  section,  the  frequency  band  and  spacing  of  vicinity 
rays  is  such  that  the  paraxial  Fresnel  volume  is  located  between  the  central 
ray  and  the  vicinity  ray  except  near  the  x-caustic  and  y-caustics  as  defined 
in  Chapman  and  Drummond  (1982).  In  the  vicinity  ray  tracing  system,  the 
x-caustic  corresponds  to  —  0  and  the  y-caustic  —  0. 

The  physical  meaning  of  is  the  distance  from  the  central  ray  to  the 
vicinity  ray,  and  is  determined  fully  from  equations  (31)  and  (32).  The 
variation  of  <7,  describes  the  change  in  amplitude,  and  the  variation  of  i], 
describes  the  geometry  of  the  wavefront.  These  properties  of  q,  and  r/,  ran  be 
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applied  to  many  related  problems  such  as  two  point  ray  tracing,  calculation 
of  the  travel  time  of  a  vicinity  ray,  amplitude  inversion,  and  correction  of 
qi  and  p;  at  a  discontinuity.  The  vicinity  ray  tracing  system  in  equations 
(31)  and  (32)  is  obtained  without  using  any  paraxial  approximations.  Thus, 
the  system  should  give  more  accurate  predictions  of  the  wavefield  in  the 
neighborhood  of  a  central  ray  when  the  medium  has  strong  velocity  gradients. 

As  shown  in  equations  (31)  and  (32),  the  vicinity  ray  tracing  system 
is  described  by  4  equations  in  3-1);  by  contrast,  the  standard  dynamic  ray 
tracing  system  requires  8  equations.  The  reason  that  the  vicinity  ray  tracing 
system  requires  fewer  equations  is  that  only  ray  centered  coordinates  are  used 
instead  of  a  combination  of  ray  centered  coordinates  and  ray  coordinates. 


2.2.3  Computation  of  travel  time  near  a  central  ray 

The  computation  of  the  travel  time  to  a  receiver  near  a  central  ray  is  just 
as  simple  in  the  vicinity  ray  tracing  system  as  in  the  standard  method  of 
dynamic  ray  tracing  using  the  paraxial  approximation.  Figure  3  illustrates 
the  calculation  of  the  travel  time,  r(a,ni,n2)  at  point  N(s,7i1,n2).  The  de¬ 
termination  of  s  and  rc,  for  a  specified  point  N  in  ray  centered  coordinates 
is  important  to  obtaining  accurate  estimates  of  travel  time  and  amplitude  of 
the  vicinity  ray  with  respect  to  a  central  ray.  The  rough  approximations  con¬ 
tained  in  the  standard  paraxial  technique  may  produce  spurious  oscillations 
in  the  superposition  of  Gaussian  beams  (e.g.,  Madariaga,  1984)  and  break 
down  if  the  central  ray  is  far  from  the  receiver. 

Here,  we  describe  an  improved  method  for  the  determination  of  a  spec¬ 
ified  point  in  ray  centered  coordinates.  We  begin  by  writing  travel  time 
field,  r(a,ni,n2)  of  the  specified  point  (e.g.,  receiver),  N  in  the  ray  centered 
coordinates  as 


r(s,ni,n2)  -  t(s)  +  £ 


where 


t j  -  l  for  qn>  x  q,  >  0  :  convex  wave  front. 

Ei  0  for  qn<  x  qi  0  :  planar  wave  front. 
c,  —  1  for  qrit  x  i/i  <  0  :  concave  wave  front. 

c,  corresponds  to  the  sign  of  M  in  the  standard  dynamic  ray  tracing  system 
along  the  (Cervcny  and  Pscncik,  1983).  At,  is  the  travel  time  difference 
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along  the  t,  direction  between  the  points  N  and  S  (figure  3).  The  travel  time 
difference  Ar(s,qni ,  9^)  between  S  and  N,  is  obtained  from 


where  An;  is  the  distance  between  C  and  N.  The  distance  An,  is  simply 
calculated  as  shown  in  figure  3,  and  is  given  by 

An;  -  yfRJ  1  n,2  1U  (38) 

where  /?;  is  the  radius  of  curvature  of  the  wavefront,  which  is  perpendicular 
to  e;.  At;  is  obtained  by  substituting  equation  (38)  into  (37), 

An  ...**  £*'*  '<•  (39) 

u(.s)  u(.s) 

To  facilitate  comparison  with  the  standard  Gaussian  beam  method,  equation 
(39)  can  be  expanded;  for  ^  -C  1, 


R,  .  n2  n2 

(,0) 

The  factor  in  equation  (40)  corresponds  to  real  M  in  the  expression  for 
the  standard  Gaussian  beam  method,  where  n,  is  small.  The  travel  time 
r(s,ni,n2)  can  be  rewritten  by  substituting  equation  (39)  into  (36), 

*  J R 2  f  n2  Ri 

r(.s,nt,n2)  -  r(j)  j  ^  r, - - -  (41) 

FI  v(*) 


Note  in  equation  (41)  that  the  travel  time  of  a  specified  point  (.s,rii,n2)  is 
obtained  without  using  any  paraxial  approximations.  The  travel  time  of  a 
specified  point  r(s,n],n2)  is  easily  and  accurately  calculated  with  respect 
to  the  travel  time  of  central  ray  r(j)  since  the  radius  of  curvature  of  the 
wavefront,  /{;,  is  a  function  of  9;  and  r/,. 

Let  9/;  denote  the  normal  distance  to  the  central  ray  from  point  1),  where 
B  is  the  intersection  of  the  wavefront  of  the  central  ray  and  vicinity  ray 
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(figure  3).  The  curvature  (Ki)  ami  the  radius  of  curvature  of  the  wavefront 
(Ri)  can  be  expressed  as  follows  (Cervcny,  1981). 

«,  vpt 


Substituting  equation  (29)  into  (42)  produces, 

„  n'x  I'*  <l’i  <?. 

Ri 

vp,  sui  7,  sm  7,  tan  T), 

Let  7,  denote  the  distance  from  S  to  B  along  the  wavefront.  The  wavefront 
coordinates  (3,71,72)  are  defined  from  projections  of  the  unit  vectors  /;  to 
the  plane  normal  to  the  central  ray  that  become  the  unit  vectors  of  the  ray 
centered  coordinates  C\  and  e2  (figure  3).  The  relation  between  7^  anf  qt  can 
be  represented  using  R,: 

R,  -  Si  (44) 

tan  rji  rji 

Equation  (44)  determines  the  Jacobian  J  between  the  ray  centered  coordi¬ 
nates  and  the  wavefront  coordinates: 


7172  -  JqiQ2 


(45) 


where  J  -  - — UX-32 .  Equation  (44)  shows  that  the  curvature  of  a  wavefront 

Ki  or  the  radius  of  curvature  of  the  wavefront  Rt  can  be  written  as  a  simple 
function  of  71  and  71-  When  n<  -  qi ,  equation  (39)  can  be  rewritten  by  using 
equation  (44), 


<7. 


v{s)tan  7, 


(^1  I  tan2T]i  l) 


(46) 


By  substituting  equation  (46)  into  (36),  the  travel  time  of  the  vicinity  ray 
Tisx(hx(ti)  can  be  expressed  as: 


"(3>  7«)  -  t(s)  I  TIT1 - 

v(s)tan  7, 


( 0  4  tan 2i),  -  1) 


(47) 


Equation  (47)  shows  that  At  also  can  be  calculated  by  just  using  7,  and  7 ; 
without  calculating  Ki  or 

r(3,n,)  and  r(s,qt)  in  equations  (41)  and  (47)  can  be  described  in  terms 
of  a  known  point  9  -  .i„  along  the  central  ray.  The  quantity  r(s)  can  be 


15 


can  be  expanded  with  respect  to  r(a„)  by  using  a  Taylor  expansion  about  a0. 
Terms  higher  titan  second  order  are  negligible  and  will  be  neglected. 

r(5)  ~  r(5o)  I  {^~  I.  ,.  (*  *o)  I  ^ fa*-*  i«  •»  (s  s<> )Z  1  -  (48) 

It  is  easy  to  see  that 


<M-S) 


and  that 


^1...  ^  («) 

da2  v 2(a0) 

The  travel  time  of  the  central  ray  r(.s)  can  be  rewritten  by  substituting 
equation  (49)  into  (48). 

r(.s)  t(.s„)  \  1  (a  ,sn )  an)2  (50) 

t’K)  2  v2(a„) 

Combining  the  expressions  (41),  (47),  and  (5U),  the  travel  time  r(.s,n.)  and 

r(.s,ql)  is  approximated  by 

/  \  /  ,  1  ,  \  1  (,si>)  /  ,2 

r(s,nt)  --  r(.s„  1  ——  (a  a„)  - (.s  .s0 

u(M  2  v (•■»..) 

*  //.) 

— 

/  .  /  ,  ,  1  /  \  1  V  1‘  {so)  ,  .2 

T(«.9.  r(-so)  '  ~r~\  's  n  ~Ti — \  ('s  ,s") 

v(aQ)  2  v2(a0) 

2 

*  Y.  Tir-  (V1  1  lanhi'  ')  (S') 

it  v(a)tan  q, 


r(s,ql)  denotes  the  travel  time  of  the  vicinity  ray,  and  q ,  is  calculated  by 
equations  (31)  and  (32),  while  r(s,u,)  indicates  the  travel  time  of  a  specified 
point  N,  such  as  a  receiver  point,  in  the  ray  centered  coordinates.  Note  that 
although  a  Taylor  expansion  has  been  used,  it  is  a  laylor  expansion  along 
the  direction  of  the  central  ray  rather  than  along  a  direction  perpendicular 
to  the  central  ray.  The  standard  Gaussian  beam  and  paraxial  ray  methods 
make  a  Taylor  expansion  of  travel  time  in  the  direction  perpendicular  to  the 
central  ray  as  well.  In  equation  (51)  it  is  usually  possible  to  select  ,s„  $, 

avoiding  the  Taylor  expansion  along  the  rrntra)  rav. 
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2.2.4  Fresnel  beamwidths 


The  bearnwidth  h\  is  defined  as  the  distance  from  the  central  ray  to  the 

paraxial  Fresnel  volume  along  the  e,  direction.  The  Fresnel  volume  encloses 

all  virtual  rays  in  the  neighborhood  of  the  central  ray  such  that  the  travel 
time  of  any  virtual  ray  differs  from  the  travel  time  of  the  central  ray  by  one- 
half  period  (Stone,  1963;  Kravtsov  and  Orlov,  1980;  Marcuse,  1982;  Cerveny, 
1988).  Given  a  frequency,  the  Fresnel  volume  is  estimated  from  the  locus  of  a 
vicinity  ray  and  a  central  ray.  The  travel  time  difference  between  the  central 
ray  and  the  paraxial  Fresnel  volume  is  given  as 

r(.s,  /•,,<))  r(s,  0,0)  7 

r(.s,0,  F2)  ~  r(.s,0,0)  7  (52) 

where  7  is  the  half  period.  Equation  (52)  states  that  a  point  on  the  Fresnel 
volume  ( .s ,  ,  0 )  or  (s,0,  b\)  has  a  half  period  time  difference  with  respect 

to  the  wavefront  that  passes  through  the  point  (5, 0,0)  on  the  central  ray. 
Formulae  for  the  beamwidths  l\  are  then  obtained  by  substituting  (51)  and 
(44)  into  (52), 

Cl  -  y/yi?  -I  27v, ■-  /W  f  27 Vi  W, 

>\  ~  . /i’I?  I  I  27 V2R7  (53) 

For  high  frequency  (small  7),  the  bearnwidth  given  by  (53)  is  approximately 
proportional  to  li\  Vx  \J ■  Equation  (53)  is  the  same  as  the  classical 
definition  of  Fresnel’s  half  period  zones  (e.g.,  Jenkin’s  and  White,  1937). 

2.3  The  synthesis  of  seismograms 

The  zeroth  order  high  frequency  asymptotic  solution  to  the  wave  equation  in 
generalized  coordinates  was  given  in  equation  (14).  bet  us  consider  a  point 
N,  located  close  to  the  central  ray,  specified  by  the  ray  centered  coordinates, 
(.9,  n( ,  n2),  that  is,  S-  (.9,  n | ,  n2).  I’he  zeroth  order  asymptotic  solution  to  the 
reduced  wave  equation  in  the  ray  centered  coordinates  can  be  expressed  in 
the  form 

,I(.S> (  .)V.M  (54) 
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The  amplitude  function  A  is  real  and  can  he  determined  hy  applying  the 
conservation  of  energy  (lux  and  a  normalization  condition  (Beiscr,  1969,  p. 
156;  Casiorowicz,  1974,  p.  45).  r  is  the  travel  time  of  the  central  ray.  k  is 
the  value  of  the  KM  AH  index  whose  value  is  increased  hy  one  whenever  the 
sign  of  qi  changes  along  the  ray.  The  KM  A  II  index  represents  the  *  phase 
shift  whenever  the  ray  touches  a  caustic  surface  (Chapman  and  Drummond. 
1982). 


2.3.1  Source-time  functions 

For  a  source-time  function  f(t)  specified  as  the  real  part  of  an  analytic  func¬ 
tion  y(t),  the  wave  field  is  given  hy  evaluating  a  convolution: 

u(S,f)  -  Hc\g{S,t)  *  y(f)l  (55) 

Rewriting  equation  (55)  as 

fir  ra o 

u(S,t)  ~  /  y(S,u>)i/(u>)  e  dej  (56) 

7T  J  oc 

and  substituting  ecpiation  (54)  for  y(5,u>)  gives 

u(S,t)  -  lie  (  (  i)k  xgn(u))y(u>)(  lu(t  r)dw  lie\y{t  r)(  i)*| 

7T  J  oc  7T 

(57) 

y(t)  is  the  analytic  time  series  represented  by 

y(l)  /(<)  ih(t)  (58) 

where  f(t)  and  h(t)  are  Hilbert  transform  pairs  (Bracewell,  1978).  'I  he  an¬ 
alytic  function  corresponding  to  any  realistic  source- time  function  can  he 
constructed  by  choosing  y(t)  to  he  a  generalized  delta  function  and  convolv¬ 
ing  that  function  with  y(t).  Some  possible  forms  for  y{t)  are 

1.  a  delta  function  (Chapman, 1 977;  Chapman  and  Drummond,  1982) 

y(t)  -  b(l.)  ij-t  (59) 


2.  a  Gaussian  wavelet 


1  («p 

—  r  v  i ' 
7T  *f 


(60) 
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3.  or  a  resonance  function  (Madariaga  and  Papadimitriou,  1985) 


tA2  i  (A t)2  "  V  l  (At)2 

7  is  the  prevailing  frequency  of  the  Gaussian  wavelet  and  A  is  the  sampling 
interval  or  the  half  period  for  the  discrete  time  series  y(t). 

The  Gaussian  wavelet  (60)  is  useful  for  simulating  a  narrow  hand  source, 
while  (61 )  is  useful  for  simulating  broad  band  responses.  Equations  (59),  (60) 
and  (61)  can  be  constructed  to  be  a  generalized  delta  function  by  requiring 

(62) 


f  He  y(t)dt 

J  no 


2.3.2  Beamwidths 

For  a  point  source,  it  is  natural  to  assume  that  the  amplitude  distribution 
within  a  beam  is  Gaussian.  The  amplitude  function  A  at  a  specified  point  N 
will  be  described  as  a  generalized  Gaussian  function  of  the  form 

MS)  crj-p|  (£-)’  (^)’l  (63) 

where  F,  are  the  half  widths  of  the  paraxial  volume.  With  this  amplitude 
distribution,  the  energy  (probability  of  finding  a  ray)  along  the  the  paraxial 
Fresnel  volume  of  half  width  /•’,  is  proportional  to  ^  and  its  amplitude  is 
proportional  to  £  with  respect  to  the  central  ray.  If  the  expression  for  the 
amplitude  function  is  viewed  as  a  generalized  delta  function,  C  can  be  chosen 
from  a  normalization  condition  in  space  and  time.  Generally,  /*’,  is  chosen 
to  be  equal  to  the  half-width  of  the  paraxial  volume  as  defined  in  equation 
(53).  We  assume  that  the  contribution  of  a  beam  to  the  receiver  is  zero 
when  the  distance  between  the  ray  and  receiver  is  larger  than  qt.  Without 
loss  of  generality,  this  condition  eliminates  a  singularity  near  both  x  and  y 
caustics.  Near  the  x  and  y  caustic,  the  half-width  of  the  paraxial  Fresnel 
volume  F,  may  become  larger  than  qt.  In  the  these  regions  the  beamwidth  is 
taken  to  be  equal  to  q,  rather  than  Ft  and  its  amplitude  distribution  will  be  a 
generalized  rectangular  function  instead  of  a  generalized  Gaussian  function. 
These  modifications  near  x  and  y  caustics  do  not  violate  energy  conservation 
and  the  normalization  condition. 
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2.3.3  Energy  conservation  and  normalization  of  beams 

Using  equations  (54)  and  (57)  the  complex  displacement  u  of  a  beam  specified 
at  a  point  N  in  the  ray  centered  coordinates  is 

«(S,l)-C«|.|-(£),-(£),|  (  .)%((■  r),  (64) 

Cl  r  2 

where  amplitude  factor  C  is  obtained  by  using  the  law  of  conservation  of 
energy  flux  and  a  normalization  condition.  This  approach  for  determining  C 
differs  from  the  approach  followed  in  the  standard  Gaussian  beam  technique, 
where  C  is  obtained  by  evaluating  the  superposition  integral  by  stationary 
phase  and  requiring  the  result  to  reproduce  ray  theory  in  regions  where  it 
is  valid.  In  contrast  to  the  standard  Gaussian  beam  method,  beams  are 
interpreted  as  the  probability  of  finding  a  ray  at  a  given  point  and  time. 
This  probability  distribution  is  assumed  to  be  a  Gaussian  distribution  whose 
unit  area  is  always  1  with  respect  to  n,  and  t.  This  constraint  guarantees 
that  the  energy  of  a  beam  (ray  tube)  is  conserved  with  respect  to  the  space 
and  time,  and  that  the  determinant  of  the  propagator  matrix  in  dynamic 
ray  tracing  system  is  constant  along  the  ray  ( biouville’s  theorem)  (Oerveny 
and  Psencik,  1983;  Kim,  1985;  Klimes,  1988).  The  wave  function  u(S,t)  then 
describes  the  probability  of  finding  a  ray  with  a  statistical  state,  which  is 
characterized  by  u.  Because  the  total  energy  in  a  beam  is  conserved  along 
the  wavefront,  it  is  necessary  to  transform  from  the  ray  centered  coordinates 
to  the  wavefront  coordinates.  The  Jacobian  between  the  ray  centered  and 
the  wavefront  coordinates  is  given  in  equation  (45).  The  conservation  law  of 
energy  flux  and  the  normalization  condition  imply  that  the  probability  (P) 
of  finding  a  ray  within  a  given  space  is 

P{S,t)  -<  u  I  u  >,,0,0'  /  /  di/i  d(f7  J  ’•  r  u‘(Sj)  1  (65 ) 


where  the  symbol  *  denotes  the  complex  conjugate  The  constant  (’  in  equa 
tion  (64)  is  determined  by  solving  equation  (65)  and  by  lonsidcring  the  nor 
maiization  factor  1)  for  a  radiation  pattern  at  the  source.  The  integral  of 
equation  (65)  yields 


P(S,t) 


a 


,2  /  j  /••>*  J 


2 


1 

D 2 


(66) 
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The  constant  [)  depends  on  the  take-off  angle  ft  and  azimuth  (f)  (Aki  and 
Richards,  1980,  p.  82;  Klimes,  1984b).  The  constant  C  is  then  obtained  as 


C 


y/2D 


VJrTf]  /‘2 

A(S)  can  be  rewritten  as  follows  by  substituting  expression  for  C  into  (63) 


MS)  - 


\/2 1) 


I  rn'\2  /  7l2  n2| 

WexP\  (y)  ■■  (y)  1 


(67) 


s/J  7T  Ty  f\ 

The  calculation  of  the  travel  time  r(.s,n,)  in  (54)  is  given  in  equations  (51). 
The  final  wavefield  of  a  ray  at  a  specified  point  is  is  then  obtained  as 


u(S,L) 


\frU) 

s/.i  7T3  /’,  /' 


tau)  ( --*)* 


(68) 


Equation  (68)  can  be  modified  to  represent  source-time  functions  such  as 
a  generalized  rectangle  function  or  combination  of  generalized  functions  for 
a  line  source.  The  representation  of  a  line  source  depends  on  its  intensity 
distribution  with  length.  Generally,  the  expression  for  a  line  source  is  more 
complicated  than  that  for  a  point  source. 

The  form  of  equation  (68)  can  be  shown  to  be  quite  similar  to  the  ex¬ 
pression  for  a  standard  Gaussian  beam  when  paraxial  approximations  are 
substituted  for  phase  and  the  expression  of  the  half  width  of  the  paraxial 
volume  is  substituted  for  F,.  The  approach  and  concepts  used  in  deriving 
the  vicinity  ray  tracing  system,  however,  are  quite  different  from  those  used 
in  the  standard  Gaussian  beam  method.  The  number  of  equations  required 
in  this  method  is  9,  by  contrast  to  21  in  the  standard  Gaussian  beam  method. 
This  method  uses  exact  positions  of  vicinity  rays  while  the  standard  Gaus¬ 
sian  beam  method  uses  estimated  values  based  on  a  Taylor  expansion  about 
the  central  ray.  Beannwidth  in  this  method  is  the  distance  from  the  central 
ray  to  the  paraxial  Fresnel  volume  and  is  fully  determined  by  equation  (53), 
while  beamwidth  in  the  standard  Gaussian  beam  method  is  usually  chosen 
arbitrarily  and  not  given  any  physical  meaning. 


2.3.4  Superposition  of  beams 

For  the  wavefield  obtained  by  superposition  of  all  beams  we  shall  use  upper¬ 
case  U,  instead  of  lower-case  ij,  which  is  reserved  for  an  individual  beam. 
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Note  that  the  wavefield  u  corresponding  to  an  individual  beam  is  a  function 
of  vertical  take-oir  angle  and  azimuth  (6  and  </>),  which  specify  the  central 
ray  under  consideration.  Thus  we  shall  write  u(5',  t,6,<£)  instead  of  u(S,l). 
The  wavefield  U(S,t)  will  be  described  by  superposition  of  individual  beams 
(rays)  with  respect  to  6  and  <f>. 


-  /  /  u{S,t,8,(f))dbd4) 

J 6.,  J 


When  the  integrand  of  equation  (69)  is  sufficiently  smooth  for  a  given  S  and 
t,  it  can  be  discretized  as 

N  M 

V(S,t)  =  (7°) 

j-  Ok  0 

where  the  quantities  Ab}  and  are  determined  from  a  given  system  Sj 
and  (fib  (Cerveny,  1983b).  The  wavefield  is  calculated  in  (69)  or  (70)  by 
summing  up  each  ray’s  contribution  at  a  specified  point,  its  wavelet  having 
a  Gaussian  distribution  both  in  space  and  time.  As  in  the  Gaussian  beam 
method  (Cerveny,  1983),  this  method  also  does  not  require  two-point  ray 
tracing  to  compute  the  seismic  wavefield. 

Since  the  energy  of  a  beam  is  conserved  along  the  wavefront,  U(s,t)  in 
equation  (69)  can  be  rewritten  in  the  wavefront  coordinates  by  using  the 
Jacobians  in  equation  (45) 

U(S,I)  --  I ‘  \  ( -p-,)’  y('l')  (  i)"dW  (71) 

VJ  rFi/(2  J\r\  J2t2 

As  a  shown  in  equations  (45)  and  (45),  the  Jacobians  J\  and  J2  between  the 
ray  centered  and  the  wavefront  coordinates  give 

J  i  ti  |  n  j 

n,  R,n,  (72) 

where  R,  is  the  radius  of  wavefront  curvature. 


2.3.5  Superposition  in  a  homogeneous  medium 

It  is  easy  demonstrate  that  the  superposition  integral  returns  simple  ray 
theory  in  a  homogeneous  medium.  In  a  homogeneous  medium,  U(S,t)  is 
represented  simply  because  Hi  —  R?  —  5,  and  S  is  the  total  distance  from 
the  source  to  the  receiver.  The  parameters  in  a  homogeneous  medium  are 
given  as: 

R.2  -  S  R2  =  S 

a  -  Ci  +  6 

~  c2  \  <f>  (73) 

and 

da  -  db 

dfl  d<t>  (74) 

where  C|  and  c2  are  constant.  Substituting  equations  (72), (73), and  (74)  into 
(71),  then  gives 


u (S-‘>  - 

SW  /v»*\  /  \k  /-re 


SyJ 7T 


y(T)(  -i)h 


where  J  J\J2 

Equation  (75)  shows  that  the  displacement  of  U(s,t)  in  a  homogeneous 
medium  is  proportional  to  the  distance  between  source  and  receiver.  Note 
that  the  superposition  is  independent  of  the  choices  made  for  the  initial  con¬ 
ditions  for  the  spacing  of  vicinity  rays. 


3  Numerical  Example 

A  numerical  example  was  chosen  to  test,  the  superposition  of  Gaussian  beams 
defined  from  the  vicinity  ray  tracing  system.  The  model  was  not  designed  to 
be  geophysically  realistic,  but  rather  to  illustrate  the  theoretical  phenomena 
near  the  caustic.  The  velocity  of  the  medium  is  given  by 

v(z)  2.5  4  O.i  x  ^  for  z  10  Km 

v(z)  3.5  |  0.4  x  ( z  10)  for  z  '  It)  Km 

The  model  has  a  discontinuity  in  the  first  derivative  of  velocity  at  z-  10  Km 
but  velocity  itself  remains  continuous  (figure  4).  In  the  vicinity  ray  tracing 
system,  no  phase  matching  method  is  required  at  a  discontinuity  of  the  first 
derivative  of  velocity.  Figure  5  shows  the  results  of  ray  tracing  , showing  a 
triplication  in  the  range  32.4  to  48.3  Km  from  the  source.  Two  caustics  are 
located  at  x  ^32.1  Km  and  x  %4ti.3  Km. 

Figure  (i  shows  the  vertical  component  synthetic  seismograms  computed 
by  superposing  Gaussian  beams  defined  from  vicinity  rays,  called  ”VRT  seis¬ 
mograms. ”  To  calculate  these  seismograms,  64  rays  are  used  with  a  1°  incre¬ 
ment  of  take  ofT  angle.  The  source  is  assumed  to  be  an  explosive  point  source, 
and  the  initial  conditions  are  i)0  5  degrees  and  </„  0.  A  monochromatic 

pulse  of  frequency  5  llz  is  used.  The  beamwidths  are  defined  by  equation  (53) 
using  the  definition  of  the  Fresnel  volume.  VVKHJ  synthetic  seismograms  are 
computed  for  comparison  with  the  VH1  seismograms  (figure  7).  As  shown 
in  figures  7  and  8,  the  two  methods  closely  agree  with  one  another.  Am¬ 
plitude  differences  between  two  methods  are  less  than  5  %.  Diffractions  are 
shown  near  the  caustics  in  both  methods.  The  diffraction  near  the  caustic 
at  x  2s32.4  Km  decays  faster  than  near  the  caustic  at  x  ~46.4  Km  because 
the  beamwidth  (Fresnel  volume)  varies  more  slowly  at  the  former  than  at  the 
latter.  Some  differences  in  the  frequency  content  of  the  diffraction  from  the 
caustic  at  x  %32.4  km  can  be  seen.  These  differences  were  generated  by  al 
lowing  for  a  broad  frequency  spectrum  in  the  WKIU  method  by  using  a  delta 
function  source  and  then  convolving  the  result  with  a  narrow  band  source 
pulse.  The  superposition  of  vicinity  rays,  on  the  other  hand,  is  such  that 
only  the  frequencies  conta  ried  in  the  narrow  band  source  pulse  can  be  seen. 
This  is  because  the  half-width  of  the  paraxial  Fresnel  volume  (beamwidth) 
was  computed  only  for  the  center  frequency  of  the  narrow  band  source  pulse. 


4  Conclusions 


Because  paraxial  approximations  are  made,  the  standard  dynamic  ray  trac¬ 
ing  system  works  in  only  a  limited  region  near  the  central  ray  at  distances  less 
than  the  scale  length  of  the  medium.  Its  use  of  two  coordinate  systems  results 
in  a  an  increased  number  of  equations  needed  to  describe  the  dynamic  prop¬ 
erties  of  a  wavefront.  An  improved  system  of  dynamic  ray  tracing  equations 
can  be  developed  from  the  Hamiltonian  of  a  ray  and  its  canonical  equations. 
This  improved  system,  which  we  have  called  the  vicinity  ray  tracing  equa¬ 
tions,  is  specified  by  only  \  equations  in  addition  to  the  kinematic  ray  tracing 
equations.  By  contrast,  the  standard  dynamic  ray  tracing  equations  based 
on  paraxial  approximations  require  8  equations  and  their  associated  Gaus¬ 
sian  beams  require  16  equations.  Unlike  the  standard  dynamic  ray  tracing 
system,  the  vicinity  ray  tracing  system  does  not  require  the  evaluation  of 
second  spatial  derivatives  of  velocity  along  a  ray.  The  vicinity  ray  tracing 
equations  will  thus  have  advantages  in  accuracy  and  in  computer  time  over 
standard  dynamic  ray  tracing.  Since  only  first  spatial  derivatives  of  velocity 
are  used  in  the  vicinity  ray  tracing  system,  no  phase  matching  is  required  at 
discontinuities  in  velocity  gradient.  At  velocity  discontinuities,  new  initial 
conditions  on  the  vicinity  rays  are  determined  by  applying  Snell’s  Law  to  the 
vicinity  rays. 

By  applying  conservation  of  energy  flux  and  a  normalization  condition,  an 
improved  Gaussian  beam  technique  can  be  developed  by  superposing  vicinity 
rays.  In  this  superposition,  bearnwidths  are  set  equal  to  the  half-width  of 
the  paraxial  Fresnel  volume.  An  example  calculation  demonstrates  that  this 
definition  of  beamwidth  can  approximate  diffraction  effects. 
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Figure  Captions 
Figure  1 

The  ray  centered  coordinates  (s, <?),<??):  t  is  the  unit  tangent  vec¬ 
tor  of  a  centra!  ray  and  ej  and  e2  are  the  unit  normal  vectors  to  t. 
The  coordinate  s  measures  the  arclength  along  a  central  ray  from 
an  arbitrary  reference  point,  q \  and  <?j  represent  length  coordi¬ 
nates  and  form  a  two  dimensional  orthogonal  coordinate  system 
in  the  plane  normal  to  fi  at  0. 

Figure  2 

The  geometry  of  the  vicinity  ray  tracing  system:  i/,  is  the  angle 
difference  between  the  tangential  vector  i  of  a  central  ray  and  the 
tangential  vector  of  a  vicinity  rav  in  the  t  -  e,  plane,  q,  is  the 
distance  between  the  central  ray  and  the  vicinity  ray  in  the  t  —  e, 
plane  at  S. 

Figure  3 

The  ray  centered  coordinates  (s,qt)  and  the  wavefront  coordi 
nates  ( s,qt ):  The  Jacobian  J  between  two  coordinates  is  given  in 
equation  (46).  The  curvature  K,  and  the  radius  of  the  curvature 
R,  of  the  wavefront  arc  described  in  terms  of  t/,  and  qt.  q[  is  the 
normal  distance  between  B  and  C. 

Figure  4 

A  laterally  homogeneous  velocity  model.  The  gradient  of  the 
velocity  has  a  discontinuity  at  z=10  Km. 

Figure  5 

Ray  trajectories  in  the  model  shown  in  figure  4.  The  triplication 
zone  is  located  in  the  range  x  =  32. 4-48. 3  Km. 
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Figure  6 

The  vertical  component  VRT  seismograms  for  the  model.  The 
center  frequency  of  &  narrow  band  Gaussian  source  wavelet  is  10 
Hz,  and  the  receivers  are  located  at  the  surface  (z=0  Km). 

Figure  7 

The  vertical  component  WKBJ  seismograms  for  the  model  for  10 
Hz.  The  conditions  are  the  same  as  in  the  VK.T  seismograms  in 
figure  6  except  that  the  WKBJ  seismograms  were  first  synthesized 
for  a  delta  source-time  function  and  then  convolved  with  a  narrow 
band  Gaussian  wavelet. 
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Figure  6:  Synthetic  Seismograms  (5  H&,  YRT) 
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Figure  7:  Synthetic  Seismograms  (5  Hz.  WKBJ) 
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Postfach  510153 

0-3000  Hannover  51 

FEDERAL  REPUBLIC  OF  GERMANY 


Ms.  Eva  Johann isson 
Senior  Research  Officer 
National  Defense  Research  Inst. 

P.0.  Box  27322 
S-102  54  Stockholm 
SWEDEN 

Tormod  Kvaarna 
NTNF/N0RSAR 
P.0.  Box  51 

N-2007  Kje I ler,  NORWAY 

Mr.  Peter  Marshall,  Procurement 
Executive,  Ministry  of  Defense 
Blacknest,  Brimpton, 

Reading  FG7-4RS 
UNITED  KINGDOM  (3  copies) 

Dr.  Robert  North 
Geophysics  Division 
Geological  Survey  of  Canada 
1  Observatory  crescent 
Ottawa,  Ontario 
CANADA,  K 1 A  0Y3 

Dr.  Frode  R i ngda I 
NTNF/N0RSAR 
P.0.  Box  51 

N-2007  Kjeller,  NORWAY 

Dr.  Jorg  Sch I ittenhardt 

Federal  Inst,  for  Geosciences  &  Nat* I  Res. 

Postfach  510153 

D-3000  Hannover  51 

FEDERAL  REPUBLIC  OF  GERMANY 

University  of  Hawai i 
Institute  of  Geophysics 
ATTN:  Dr.  Daniel  Walker 
Honolulu,  HI  96822 
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FOREIGN  CONTRACTORS 


Or.  Ramon  Cabre,  S.J. 

Observatorio  San  Calixto 
Casi I  la  5939 
La  Paz  Bolivia 

Professor  Peter  Harjes 
Institute  for  Geophysik 
Rhur  University/Bochum 

F. O.  Box  102148,  4630  Bociium  1 
FEDERAL  REPUBLIC  OF  GERMANY 

Or.  E.  Husebye 
NTNF/NORSAR 
P.0.  Box  51 

N-2007  Kjeller,  NORWAY 

Professor  Brian  L.N.  Kennett 
Research  School  of  Earth  Sciences 
Institute  of  Advanced  Studies 

G. P.O.  Box  4 
Canberra  2601 
AUSTRALIA 

Or.  B.  Mass i non 

Societe  Radiomana 

27,  Rue  Claude  Bernard 

7,005,  Paris,  FRANCE  (2  copies) 

Dr.  Pierre  Mechler 
Societe  Radiomana 
27,  Rue  Claude  Bernard 
75005,  Paris,  FRANCE 

Or.  Svein  Mykkeltveit 

NTNF/NORSAR 

P.0.  Box  51 

N-2007  Kjeller,  NORWAY  (3  copies) 
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GOVERNMENT 


Or.  Ralph  Alewine  III 
DARPA/NMRO 

1400  Wilson  Boulevard 
Arlington,  VA  22209-2308 

Or.  Robert  Blandford 
0ARPA/NMR0 

1400  Wilson  Boulevard 
Arlington,  VA  222o9-*_>G8 

Sand i a  National  Laboratory 
ATTN:  Or.  H.  B.  Durham 
Albuquerque,  NM  87185 

Or.  Jack  Evernden 
USGS-Earthquake  Studies 
345  Middlefield  Road 
Menlo  Park,  CA  94025 

U'.S.  Geological  Survey 

ATTN:  Or.  T.  Hanks 

Nat’ I  Earthquake  Resch  Center 

345  Middlefield  Road 

Menlo  Park,  CA  94025 

Or.  James  Hannon 
Lawrence  Livermore  Mat1 I  Lab, 
P.0.  Box  808 
Livermore,  CA  94550 

Paul  Johnson 

ESS-4,  Mail  Stop  J979 

Los  Alamos  National  Laboratory 

Los  Alamos,  NM  87545 

Ms.  Ann  Kerr 
0ARPA/NMR0 

1400  Wilson  Boulevard 
Arlington,  VA  22209-2308 

Or.  Max  Koontz 
US  Dept  of  Energy/GP  5 
Forrestal  Building 
1000  Independence  Ave. 
Washington,  O.C.  20585 


Dr.  W.  H.  K.  Lee 
USGS 

Office  of  Earthquakes,  Volcanoes, 

&  Engineering 
8ranch  of  Seismology 
345  Middlefield  Rd 
Menlo  Park,  CA  94025 

Or.  William  Leith 
U.S.  Geological  Survey 
Ma i I  Stop  928 
Reston,  VA  22092 

Or.  Richard  Lewis 
Oir.  Earthquake  Engineering  and 
Geophys i cs 

U.S.  Army  Corps  of  Enqineers 
Box  631 

Vicksburg,  MS  39180 

Dr.  Robert  Masse' 

Box  25046,  Mail  Stop  967 
Denver  Federa I  Center 
Denver,  Colorado  tt0225 

Richard  Morrow 
ACOA/V I 
Room  5741 

320  21st  Street  N.W. 

Washington,  O.C.  20451 

Or.  Keith  K.  Nakanishi 

Lawrence  Livermore  National  Laboratory 

P.0.  Box  808,  L-205 

Livermore,  CA  94550  (2  copies) 

Or.  Carl  Newton 

Los  Alamos  National  Lab. 

P.0.  Box  1663 

Mall  Stop  C335,  Group  ESS-3 
Los  Alamos,  NM  87545 

Dr.  Kenneth  H.  Olsen 
Los  Alamos  Scientific  Lab. 

P.0.  Box  1663 

Mall  Stop  C335,  Group  ESS-3 
Los  Alamos,  NM  87545 

Howard  J.  Patton 
Lawrence  Livermore  National 
Laboratory 
P.0.  Box  808,  L-205 
Livermore,  CA  94550 
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Mr.  Chris  Paine 

Office  of  Senator  Kennedy 

SR  315 

United  States  Senate 
Washington,  D.C.  20510 

AFOSR/NP 

ATTN:  Colonel  Jerry  J.  Perrizo 
Bldg  410 

Bolling  AFB,  Wash  O.C.  20332-6448 
HQ  AFTAC/TT 

Attn:  Or.  Frank  F.  Pilotte 
Patrick  AFB,  Florida  32925-6001 

Mr.  Jack  Rach I i n 
USGS  -  Geology,  Rm  3  C136 
Mail  Stop  928  National  Center 
Reston,  VA  22092 

Robert  Reinke 
AFWL/NTESG 

Kirtland  AFB,  NM  87117-6008 

Or.  Byron  Ristvet 

HQ  DNA,  Nevada  Operations  Office 

Attn:  NVCG 

P.0.  Box  98539 

Las  Vegas,  NV  89193 

HQ  AFTAC/TGR 

Attn:  Or.  George  H.  Rothe 
Patrick  AFB,  Florida  32925-6001 

Oonald  L.  Springer 

Lawrence  Livermore  National  Laboratory 
P.0.  Box  808,  L-205 
Livermore,  CA  94550 

Or.  Lawrence  Turnbull 
OSWR/NEO 

Central  Intelligence  Agency 
CIA,  Room  5G48 
Washington,  O.C.  20505 

Or.  Thomas  Weaver 

Los  Alamos  National  Laboratory 

P.0.  Box  1663 

MS  C  335 

Los  Alamos,  NM  87545 
GL/SULL 

Research  Library 

Hanscom  AFB,  MA  01731-5000  (2  copies) 


Secretary  of  the  Air  Force  (SAFR0) 
Washington,  X  20330 
Office  of  the  Secretary  Defense 
OOR  4  E 

Washington,  X  20330 
HQ  DNA 

ATTN:  Technical  Library 
Washington,  X  20305 

DARPA/RMO/RETR I EVAL 
1400  W! Ison  Blvd. 

Arlington,  VA  22209 

DAFIPA/RMO/Secur i ty  Office 
1400  Wi Ison  Bl vd. 

Arlington,  VA  22209 

GL/XO 

Hanscom  AFB,  MA  01731-5000 
GL/LW 

Hanscom  AFB,  MA  01731-5000 
0ARPA/PM 

1400  Wilson  Boulevard 
Arlington,  VA  22209 

Defense  Technical 
Information  Center 
Cameron  Station 
Alexandria,  VA  22314 
(5  copies) 

Defense  Intelligence  Agency 
Directorate  for  Scientific  & 
Technical  Intelligence 
Washington,  D.C.  20301 

Defense  Nuclear  Agency/SPSS 
ATTN:  Dr.  Michael  Shore 
6801  Telegraph  Road 
Alexandria,  VA  22310 

AFTAC/CA  (ST  INFO) 

Patrick  AFB,  FL  32925-6001 

Dr.  Gregory  van  der  Vink 
Congress  of  the  United  States 
Office  of  Technology  Assessment 
Washington,  D.C.  20510 

Mr.  Alfred  Lieberman 

ACDA/V I -0A* State  Department  Building 

Room  5726 

320  -  2 ISt  Street,  NW 
Washington,  D.C.  20451 
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TACT EC 

Battel  I e  Memorial  Institute 
505  King  Avenue 

Columbus,  OH  43201  (Final  report  only) 
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